classdef cellSeparation < handle
    properties
        Het = 0;
        Hfin;
        Gain = [];
        Poles = 0;
        Zeros = 0;
        K = 0;
        alto = 0;
        d = inf;
    end
    
    methods 
        
        function this = cellSeparation(Htot)
            [zeros polos k] = zpkdata(Htot);
            this.K = k;
            zeros = cell2mat(zeros);
            polos = cell2mat(polos);
            this.Poles = polos;
            this.Zeros = zeros;
        end        
   
        function autoJoinPZ(obj)
            zeros = obj.Zeros;
            polos = obj.Poles;
            k = obj.K;
            [matriz G] = tf2sos(k*poly(zeros),poly(polos));
            obj.alto = length(matriz(:,1));
            obj.Hfin = tf(1,[2 3]);
            diferencia = obj.alto*2-length(zeros);

            for k=1:obj.alto

               num=[matriz(k,1) 0 matriz(k,3)];
               den=[matriz(k,4) matriz(k,5) matriz(k,6)];


               if den(3) == 0
                   den(3) = den(2);
                   den(2) = den(1);
                   den(1) = 0;
               end

                if num(3) == 0
                   if diferencia > 1
                       num(1) = 0;
                       num(3) = 1;
                       diferencia = diferencia-2;
                   else if diferencia == 1
                           num(1) = 0;
                           num(2) = 1;
                           diferencia = diferencia-1;
                       end
                   end
                end
                
                aux(k) = tf(num,den);%No me dejaba hacerlo directamente sobre Het no se porque :)

            end
            obj.Het = aux;
        end
        
%         function autoCascade(obj)
%             
%         end

        function autoEverything(obj,perm,w)
            obj.autoJoinPZ();
            obj.autoCascade(perm,w);
        end

        function autoCascade(obj, perm, w)

        dsel=inf;

            for k=1: length(perm)
                %disp('.')


                        [daux,gain(k,:),Het2]=evalua(obj.Het,perm(k,:),w);

                        if daux<=dsel
                            Hetsel=Het2;
                            selected=k;
                            dsel=daux;
                        end       
            end
            obj.Hfin = Hetsel;
            obj.Gain = gain(selected,:);
            obj.d = dsel;
            perm(selected,:);

        end
        
        function  [daux,gain,Haux] = evaluate(obj, Haux,distribution,w)
           
            gaindelta=1.1;

            etapas = length(distribution);
            gain = ones(etapas);
            gain = gain(1,:);
            daux=inf;






            for k=1:etapas
                for l=1:etapas 
                    if k~=l                  

                                Haux(distribution(k))=Haux(distribution(k))*indicedefidelidad;
                                Haux(distribution(l))=Haux(distribution(l))/indicedefidelidad;
                                gain(k)=gain(k)*gaindelta;
                                gain(l)=gain(l)/gaindelta;

                            while rango(Haux,etapas,w,configuracion)<=daux
                                daux=rango(Haux,etapas,w,configuracion);
                                Haux(distribution(k))=Haux(distribution(k))*indicedefidelidad;
                                Haux(distribution(l))=Haux(distribution(l))/indicedefidelidad;
                                gain(k)=gain(k)*gaindelta;
                                gain(l)=gain(l)/gaindelta;

                            end
                                Haux(distribution(k))=Haux(distribution(k))/indicedefidelidad;
                                Haux(distribution(l))=Haux(distribution(l))*indicedefidelidad;
                                gain(k)=gain(k)/gaindelta;
                                gain(l)=gain(l)*gaindelta;


                            while rango(Haux,etapas,w,configuracion)<=daux;
                                daux=rango(Haux,etapas,w,configuracion);
                                Haux(distribution(k))=Haux(distribution(k))/indicedefidelidad;
                                Haux(distribution(l))=Haux(distribution(l))*indicedefidelidad;
                                gain(k)=gain(k)/gaindelta;
                                gain(l)=gain(l)*gaindelta;
                            end 
                    end

                end

            end         

        end

        
        function d = range(obj,w,configuracion)
            H=1;
            min1=inf;
            max1=0;

            for k=1:largo

                    H=H*obj.Het(configuracion(k));
                    [MAG]=bode(H,w);


                    maxaux=max(MAG);
                    minaux=min(MAG);

                    if maxaux>max1
                        max1=maxaux;
                    end

                    if minaux<min1
                        min1=minaux;
                    end

            end

               d = max1/min1-1;
        end
        
    end
end